* simulate scenarios that vary WTP for fuel economy for 2012 and 2022 standards, which are reported in Table 7

* scenarios:

* 1) - 2) 2012 and 2022 standards with mean WTP for horsepower

* 3) - 4) 2012 and 2022 standards with mean price sensitivity

* save data set with simulation results, which are reported using simulations_report.do

* created 7-15-25


clear all
set more off
set type double
set mat 11000
log using table7_simulations, replace text

******************************
* prepare data for simulations
******************************

* macro for demographic group
global demo = "income_group old urban"

* macro for vehicle identifier
global vehicle = "make model series drive_type fuel_type body_style liter"

* macro for range of efficiency change
global T_range = 0.4

* macro for bound on price change
global price_tol = 1.25

* macro for range of CAFE shadow price change
global cafe_price_change = 10

* macros for shadow price on fuel consumption rate and conversion factors
global miles_cars = 195264
global miles_trucks = 225865
global mg_mile = 0.0085806

* carbon dioxide emissions factor (kg/gallon)
global ef = 8.61

* macro for rebound effect
global rebound = -0.1

* macros for bounds on changes in horsepower and fuel economy
global hp_range = 0.2
global mpg_range = 0.8

* macro for maximum zev price
global zevp_max = 8500

* macro for number of years between base year and simulation year
global num_years = 6

* electricity emissions rates
global elec_er = 0.12

* kwh/mile
global kwh_mile = 1/3

* macro for effect of 1 percent mpg change on fuel consumption rate
global delta_mpg = (1/25 - 1/25.25)

* macro for share of sales in zev states
global zev_share = 0.5

* load data
use "L:\Project-MaritzCX\Workspace1\Public\CAFE Ex Post Replication\Data\supply_estimation_results", clear

* compute efficiency in 2012 and 2018
gen eff = ln(mpg) - tradeoff_hp*ln(hp/weight) 
sort market_ID make model $demo $vehicle
by market_ID make model: egen wt = sum(sales)
replace wt = sales/wt
by market_ID make model: egen m_eff = sum(wt*eff)
replace m_eff = . if model==model[_n-1]
for X in num 2012 2018: gen yX = 1 if market_ID==X
for X in num 2012 2018: egen m_eff_X = mean(yX*m_eff)
gen delta_T = m_eff_2018 - m_eff_2012

* keep 2018 observations
drop eff wt m_eff* y2012 y2018
keep if market_ID==2018

* electricity emissions rates
gen elec_er = $elec_er

* create year-vehicle interactions and demographic groups
egen vehicle = group($vehicle)
egen id = group(market_ID vehicle)
egen demo = group($demo)
gen demo1 = demo==1
drop if id==. | demo==.

* log hp/wt
gen ln_hw = ln(hp/weight)

* compute ZEV subsidy (assume fraction of sales in zev states)
gen zev_req = 0.03 if market_ID==2015 | market_ID==2016
replace zev_req = 0.04 if market_ID==2017
replace zev_req = 0.045 if market_ID==2018
replace zev_req = 0 if zev_req==.
gen zev_credit = 0 if fuel_type=="Electric" & fuel_type!="Plug-in Hybrid" & market_ID>=2015
replace zev_credit = 0.01*evrange + 0.5 if fuel_type=="Electric" & market_ID>=2015
replace zev_credit = 4 if fuel_type=="Electric" & evrange>=350 & market_ID>=2015
replace zev_credit = 0.01*evrange + 0.3 if fuel_type=="Plug-in Hybrid" & market_ID>=2015
replace zev_credit = 1.1 if fuel_type=="Plug-in Hybrid" & evrange>=80 & market_ID>=2015
replace zev_credit = 0 if zev_credit==.
gen zev_subsidy = $zev_share*zev_price*(zev_credit - zev_req)

* credit multipliers
gen electric = fuel_type=="Electric"
gen phev = fuel_type=="Plug-in Hybrid"
gen pev = electric + phev
qui for X in num 1/4: gen multiplier_X = 1
qui for X in num 2 4: replace multiplier_X = 2 if electric==1
qui for X in num 2 4: replace multiplier_X = 1.6 if phev==1

* compute standards
gen a_2012 = 35.95 if truck==0
gen b_2012 = 27.95 if truck==0
gen c_2012 = 0.0005308 if truck==0
gen d_2012 = 0.006057 if truck==0
replace a_2012 = 29.82 if truck==1
replace b_2012 = 22.27 if truck==1
replace c_2012 = 0.0004546 if truck==1
replace d_2012 = 0.0149 if truck==1
gen a_2022 = 50.24 if truck==0
gen b_2022 = 37.59 if truck==0
gen c_2022 = 0.000447 if truck==0
gen d_2022 = 0.00159 if truck==0
replace a_2022 = 40.31 if truck==1
replace b_2022 = 26.02 if truck==1
replace c_2022 = 0.000499 if truck==1
replace d_2022 = 0.00436 if truck==1
qui for X in num 1 3: gen standard_X = 0.8/(min(max(c_2012*footprint + d_2012, 1/a_2012), 1/b_2012))
qui for X in num 2 4: gen standard_X = 0.8/(min(max(c_2022*footprint + d_2022, 1/a_2022), 1/b_2022))

* count number of vehicles sold by each firm			
gen one = 1
sort market_ID demo firm vehicle
by market_ID demo firm: egen count = sum(one)
by market_ID demo firm: gen vehnum = sum(one)

* predict market share and sales
sort demo vehicle
by demo: egen pr_mkt_share = sum(exp(util_new))
replace pr_mkt_share = exp(util_new)/(exp(util_used) + pr_mkt_share)
replace pr_mkt_share = 0 if pr_mkt_share==.
gen pr_sales = pr_mkt_share*mkt_size

* compute average efficiency change needed to meet new standard with technology
gen fcr_adj = 1/mpg
replace fcr_adj = 0 if fuel_type=="Electric"
replace fcr_adj = 0.5*fcr_adj if fuel_type=="Plug-in Hybrid"
egen wt = sum(pr_sales)
replace wt = pr_sales/wt
egen mean_fcr = sum(wt*fcr_adj)
qui for X in num 1/4: egen mean_fcr_X = sum(wt/standard_X)
qui for X in num 1/4: gen eff_X = (mean_fcr - mean_fcr_X)/mean_fcr
drop wt mean_fcr mean_fcr_*

* create variables for each scenario
qui for X in num 1/4: gen mc_X = mc*(1 + eff_X)^mc_fe
qui for X in num 1/4: gen mpg_X = mpg/(1 - eff_X)
qui for X in num 1/4: gen cpm_X = fuel_price/mpg_X
qui for X in num 1/4: replace mpg_X = mpg if fuel_type=="Electric"
qui for X in num 1/4: replace cpm_X = cpm if fuel_type=="Electric"
qui for X in num 1/4: gen fcr_X = 1/mpg_X
qui for X in num 1/4: gen fcr_adj_X = fcr_X
qui for X in num 1/4: replace fcr_adj_X = 0 if fuel_type=="Electric"
qui for X in num 1/4: replace fcr_adj_X = 0.5*fcr_adj_X if fuel_type=="Plug-in Hybrid"
qui for X in num 1 3: gen cafe_price_X = cafe_price
qui for X in num 2 4: gen cafe_price_X = cafe_price*3
qui for X in num 1/4: gen feebate_X = cafe_price_X*(fcr_adj_X - 1/standard_X)*multiplier_X
qui for X in num 1/4: gen price_new_X = trans_price + subsidy + mc_X - mc - zev_subsidy + feebate_X

* initialize price and create variable for new price
gen price_old = trans_price
gen price_new = trans_price
gen price_old_firm = trans_price
gen price_new_firm = trans_price

* initialize technology and performance
qui for X in num 1/4: gen T_X = delta_T + eff_X
qui for X in num 1/4: gen ln_hw_X = (ln(mpg_X) - ln(mpg) - T_X + delta_T)/tradeoff_hp + ln_hw

* initialize initial technology variables to keep track of change and cycling
gen T_old = T_1
for X in num 2/4: gen T_oldX = T_old
gen T_cycle = .
gen t_diff = .

* compute mean WTP for horsepower
gen wtp_hp = beta_hw/alpha
summ wtp_hp [aw=sales]
global mean_wtp_hp = r(mean)

* set hp/weight coefficient
for X in num 1 2: gen beta_hw_X = alpha*$mean_wtp_hp
for X in num 3 4: gen beta_hw_X = beta_hw

* compute mean price sensitivity
summ alpha [aw=sales]
global mean_alpha = r(mean)

* set price sensitivity paramater
for X in num 1 2: gen alpha_X = alpha
for X in num 3 4: gen alpha_X = $mean_alpha

* new utility
qui for X in num 1/4: gen util_new_X = util_new + alpha_X*(price_new_X - subsidy - trans_price) + beta_cpm*(cpm_X - cpm) + beta_hw_X*(ln_hw_X - ln_hw)

* compute share of used in total sales
gen used_sh = used_sales/mkt_size

* utility of outside opion
qui for X in num 1/4: gen util_used_X = .

* predict market share and sales
sort demo vehicle
qui for X in num 1/4: by demo: egen pr_mkt_share_X = sum(exp(util_new_X))
qui for X in num 1/4: replace util_used_X = ln(used_sh*pr_mkt_share_X/(1 - used_sh))
qui for X in num 1/4: replace pr_mkt_share_X = exp(util_new_X)/(exp(util_used) + pr_mkt_share_X)
qui for X in num 1/4: replace pr_mkt_share_X = 0 if pr_mkt_share_X==.
qui for X in num 1/4: gen pr_sales_X = pr_mkt_share_X*mkt_size
sort vehicle demo
qui for X in num 1/4: by vehicle: egen vpr_sales_X = sum(pr_sales_X)

* compute markup (difference between price, marginal costs, and feebate)
qui for X in num 1/4: gen markup_X = price_new_X - mc_X - feebate_X + zev_subsidy

* create variables for summation and revenue change
gen summation_h = .
gen summation_m = .
gen revenue_change = .

* dummy if firm only has electrics
sort firm vehicle demo
by firm: egen minelectric = min(electric)

* variables for imposing constraints on fuel economy and horsepower changes
gen constrained_high_mpg = .
gen constrained_low_mpg = .
gen constrained_high_hp = .
gen constrained_low_hp = .
foreach kind in h1_high h2_high f1_high f2_high h1_low h2_low f1_low f2_low {
	gen `kind' = .
}

* variables for own and cross attribute derivatives
gen own_deriv_h = .
gen own_deriv_m = .
gen cross_deriv_h = .
gen cross_deriv_m = .

* number of iterations
gen num_iter_f = 0
gen num_iter_z = 0
gen num_iter_t = 0

******************************
* simulations
******************************

************
* cycle through scenarios
************

* number of firms
egen firm_num = group(firm)

* macros
summ demo
global maxdemo = r(max)
summ firm_num
global max = r(max)
global tolerance = 10^-2
global price_change = 1
global firm_price_change = 1
global firm_t_change = 1
global cafe = 1
global zev = 1
global max_iter = 20

forvalues scenario = 1/4 {

************
* First loop: iterate until market-wide technology changes are small
************

	while $price_change>$tolerance {

************	
* Second loop: iterate until CAFE standard is achieved
************
		while $cafe>=$tolerance {
		    
************
* Third loop: iterate over firms and regions
************	

			forvalues num = 1/$max {
			    
************
* Loop 4a: iterate until firm's technology choices converge
************				
				
				while $firm_t_change>=$tolerance {

* increase number of iterations for technology choice
					qui replace num_iter_t = num_iter_t + 1				
				
* determine number of vehicles sold by company
					summ count if firm_num==`num'
					global numveh = r(max)					
		
* compute own derivative term
					qui replace own_deriv_h = beta_hw_`scenario'*pr_mkt_share_`scenario'*(1 - pr_mkt_share_`scenario')
					qui replace own_deriv_m = beta_cpm*fuel_price*pr_mkt_share_`scenario'*(1 - pr_mkt_share_`scenario')	
			
* compute cross derivative term
					qui replace cross_deriv_h = -beta_hw_`scenario'*pr_mkt_share_`scenario'
					qui replace cross_deriv_m = -beta_cpm*fuel_price*pr_mkt_share_`scenario'
	
* create vectors of own and cross price derivatives
					sort firm vehicle demo
					mkmat own_deriv_h if firm_num==`num'
					mkmat cross_deriv_h if firm_num==`num'		
					mkmat own_deriv_m if firm_num==`num'
					mkmat cross_deriv_m if firm_num==`num'	
					mkmat markup_`scenario' if firm_num==`num' & demo==1	

* create vector for marginal costs and feebate	
					mkmat mc_`scenario' if firm_num==`num' & demo==1
					mkmat feebate_`scenario' if firm_num==`num' & demo==1		
					mkmat zev_subsidy if firm_num==`num' & demo==1		

* create matrix with market shares	
					sort firm vehnum demo
					mkmat pr_mkt_share_`scenario' if firm_num==`num'
					matrix define S = J($numveh, $maxdemo, 0)
					forvalues j = 1/$numveh {
						forvalues k = 1/$maxdemo {
							matrix define S[`j',`k'] = pr_mkt_share_`scenario'[`k' + (`j' - 1)*$maxdemo,1]
						}
					}			

* create vector for market size
					mkmat mkt_size if firm_num==`num' & vehnum==1
	
* create matrices of derivatives (horsepower matrix used for technology choice and both matrices used for fuel consumption rate choice)
					matrix define deriv_h = J(1,$maxdemo,0)
					matrix define D_h = J($numveh, $numveh,0)
					matrix define deriv_m = J(1,$maxdemo,0)
					matrix define D_m = J($numveh, $numveh,0)
					forvalues j = 1/$numveh {
						forvalues k = 1/$numveh {
							if `j'==`k' {
								forvalues l = 1/$maxdemo {
									matrix define deriv_h[1,`l'] = own_deriv_h[(`j' - 1)*$maxdemo + `l',1]
									matrix define deriv_m[1,`l'] = own_deriv_m[(`j' - 1)*$maxdemo + `l',1]	
								}
							}
							if `j'!=`k' {
								forvalues l = 1/$maxdemo {
									matrix define deriv_h[1,`l'] = cross_deriv_h[(`k' - 1)*$maxdemo + `l',1]*pr_mkt_share_`scenario'[(`j' - 1)*$maxdemo + `l',1]
									matrix define deriv_m[1,`l'] = cross_deriv_m[(`k' - 1)*$maxdemo + `l',1]*pr_mkt_share_`scenario'[(`j' - 1)*$maxdemo + `l',1]
								}			
							}
							matrix define D_h[`j',`k'] = deriv_h*mkt_size
							matrix define D_m[`j',`k'] = deriv_m*mkt_size
						}
					}
* multiply matrix by markup
					matrix define summation_h = D_h*markup_`scenario'
					matrix define summation_m = D_m*markup_`scenario'

* fill in summation variable
					sort firm vehnum demo
					forvalues num_m = 1/$numveh {
						qui replace summation_h = summation_h[`num_m',1] if firm_num==`num' & vehnum==`num_m'
						qui replace summation_m = summation_m[`num_m',1] if firm_num==`num' & vehnum==`num_m'
					}	
				
* compute operating profits	
					sort make model vehicle demo
					by make model: egen sumprofits = sum((fcr_adj_`scenario'*(cafe_price_`scenario'*vpr_sales_`scenario' - summation_m) - vpr_sales_`scenario'*mc_fe*mc_`scenario')*demo1)

* compute technology choice (accounting for length of steady state), imposing constraints
					qui replace T_`scenario' = (sumprofits*$num_years)/(2*gamma_m) if firm_num==`num' 
					qui replace T_`scenario' = delta_T*(1 + $T_range) if T_`scenario'>delta_T*(1 + $T_range)
					qui replace T_`scenario' = delta_T*(1 - $T_range) if T_`scenario'<delta_T*(1 - $T_range)
					qui replace T_`scenario' = 0 if electric==1
					qui replace T_`scenario' = 0.5*T_`scenario' if phev==1

* compute new marginal costs
					qui replace mc_`scenario' = mc*exp(mc_fe*(T_`scenario' - delta_T)) if firm_num==`num'
					drop sumprofits

* compute fuel consumption rate from first order condition
					qui replace fcr_`scenario' = summation_h/(tradeoff_hp*(summation_m - cafe_price_`scenario'*vpr_sales_`scenario')) if firm_num==`num' & electric==0

* compute horsepower using tradeoff and technology constraint 
					qui replace ln_hw_`scenario' = ln(hp/weight) - (1/tradeoff_hp)*(ln(fcr_`scenario'*mpg) + T_`scenario' - delta_T) if firm_num==`num' & electric==0

* determine if attribute choices are constrained relative to 2011 values
					qui replace constrained_high_hp = ln_hw_`scenario'>ln(hp/weight) + $hp_range if firm_num==`num'
					qui replace constrained_low_mpg = ln(fcr_`scenario'*mpg)>$mpg_range if firm_num==`num'
					qui replace constrained_high_mpg = ln(fcr_`scenario'*mpg)<-$mpg_range if firm_num==`num'
					qui replace constrained_low_hp = ln_hw_`scenario'<ln(hp/weight) - $hp_range if firm_num==`num'

* determine levels of horsepower and fuel consumption rate if each constraint binds						
					qui replace h1_high = ln(hp/weight) + $hp_range if firm_num==`num'
					qui replace h2_high = ln(hp/weight) - (1/tradeoff_hp)*(ln(fcr_`scenario'*mpg) + T_`scenario' - delta_T) if firm_num==`num'	
					qui replace f1_high = (1/mpg)*exp((-tradeoff_hp)*(ln_hw_`scenario' - ln(hp/weight) - T_`scenario' + delta_T)) if firm_num==`num'
					qui replace f2_high = (1/mpg)*exp(-$mpg_range) if firm_num==`num'
					qui replace h1_low = ln(hp/weight) - (1/tradeoff_hp)*(ln(fcr_`scenario'*mpg) + T_`scenario' - delta_T) if firm_num==`num'
					qui replace h2_low = ln(hp/weight) - $hp_range if firm_num==`num'
					qui replace f1_low = (1/mpg)*exp(-$mpg_range) if firm_num==`num'
					qui replace f2_low = (1/mpg)*exp((-tradeoff_hp)*(ln_hw_`scenario' - ln(hp/weight) - T_`scenario' + delta_T)) if firm_num==`num'
			
***						
* impose constraints on change in fuel consumption rate and horsepower
* high horsepower binds
					qui replace ln_hw_`scenario' = ln(hp/weight) + $hp_range if firm_num==`num' & constrained_high_hp==1 & electric==0					
					qui replace fcr_`scenario' = (1/mpg)*exp((-tradeoff_hp)*(ln_hw_`scenario' - ln(hp/weight)) - T_`scenario' + delta_T) if firm_num==`num' & constrained_high_hp==1 & electric==0

* low fuel economy binds
					qui replace fcr_`scenario' = (1/mpg)*exp(-$mpg_range) if firm_num==`num' & constrained_low_mpg==1 & electric==0						
					qui replace ln_hw_`scenario' = ln(hp/weight) - (1/tradeoff_hp)*(ln(fcr_`scenario'*mpg) + T_`scenario' - delta_T) if firm_num==`num' & constrained_low_mpg==1 & electric==0

* high fuel economy binds
					qui replace fcr_`scenario' = (1/mpg)*exp(-$mpg_range) if firm_num==`num' & constrained_high_mpg==1 & electric==0
					qui replace ln_hw_`scenario' = ln(hp/weight) - (1/tradeoff_hp)*(ln(fcr_`scenario'*mpg) + T_`scenario' - delta_T) if firm_num==`num' & constrained_high_mpg==1 & electric==0		
					
* low horsepower binds
					qui replace ln_hw_`scenario' = ln(hp/weight) - $hp_range if firm_num==`num' & constrained_low_hp==1 & electric==0
					qui replace fcr_`scenario' = (1/mpg)*exp((-tradeoff_hp)*(ln_hw_`scenario' - ln(hp/weight)) - T_`scenario' + delta_T) if firm_num==`num' & constrained_low_hp==1 & electric==0

* high horsepower and low fuel economy bind
					qui replace ln_hw_`scenario' = min(h1_high, h2_high) if firm_num==`num' & constrained_high_hp==1 & constrained_low_mpg==1 & electric==0
					qui replace fcr_`scenario' = (1/mpg)*exp((-tradeoff_hp)*(ln_hw_`scenario' - ln(hp/weight)) - T_`scenario' + delta_T) if firm_num==`num' & constrained_high_hp==1 & constrained_low_mpg==1 & electric==0
						
* low horsepower and high fuel economy bind
					qui replace ln_hw_`scenario' = max(h1_low, h2_low) if firm_num==`num' & constrained_low_hp==1 & constrained_high_mpg==1 & electric==0
					qui replace fcr_`scenario' = (1/mpg)*exp((-tradeoff_hp)*(ln_hw_`scenario' - ln(hp/weight)) - T_`scenario' + delta_T) if firm_num==`num' & constrained_low_hp==1 & constrained_high_mpg==1 & electric==0						

* check that fuel consumption rate and horsepower don't change for electrics	
					qui replace fcr_`scenario' = (1/mpg) if electric==1
					qui replace ln_hw_`scenario' = ln_hw if electric==1					
***

* compute new fuel economy, cost-per-mile, performance, feebate
					qui replace mpg_`scenario' = 1/fcr_`scenario' if firm_num==`num' & electric==0
					qui replace cpm_`scenario' = fuel_price*fcr_`scenario' if firm_num==`num'
					qui replace fcr_adj_`scenario' = fcr_`scenario' if firm_num==`num'
					qui replace fcr_adj_`scenario' = 0 if electric==1 & firm_num==`num'
					qui replace fcr_adj_`scenario' = 0.5*fcr_adj_`scenario' if phev==1 & firm_num==`num'
					qui replace feebate_`scenario' = cafe_price_`scenario'*(fcr_adj_`scenario' - 1/standard_`scenario') if firm_num==`num'

* update feebate and markup
					qui replace markup_`scenario' = price_new_firm - mc_`scenario' - feebate_`scenario' + zev_subsidy if firm_num==`num'
					
* update utility and market shares
					drop util_new_`scenario' pr_mkt_share_`scenario' pr_sales_`scenario' vpr_sales_`scenario'
					gen util_new_`scenario' = util_new + alpha_`scenario'*(price_new_firm - subsidy - trans_price) + beta_cpm*(cpm_`scenario' - cpm) + beta_hw_`scenario'*(ln_hw_`scenario' - ln_hw)
					sort demo vehicle
					by demo: egen pr_mkt_share_`scenario' = sum(exp(util_new_`scenario'))
					qui replace util_used_`scenario' = ln(used_sh*pr_mkt_share_`scenario'/(1 - used_sh))
					qui replace pr_mkt_share_`scenario' = exp(util_new_`scenario')/(exp(util_used_`scenario') + pr_mkt_share_`scenario')
					qui replace pr_mkt_share_`scenario' = 0 if pr_mkt_share_`scenario'==.
					gen pr_sales_`scenario' = pr_mkt_share_`scenario'*mkt_size
					sort vehicle demo
					by vehicle: egen vpr_sales_`scenario' = sum(pr_sales_`scenario')

* compute difference between current and previous technology
					egen firm_wt = sum(pr_sales_`scenario') if firm_num==`num'
					qui replace firm_wt = pr_sales_`scenario'/firm_wt if firm_num==`num'
					egen firm_fcr = sum(firm_wt*T_`scenario') if firm_num==`num'
					egen firm_fcr_old = sum(firm_wt*T_old) if firm_num==`num'
	
* compute change in fuel consumption rate
					egen m_t_diff = max(abs((T_`scenario' - T_old)/T_old)) if firm_num==`num' & electric==0
					qui replace m_t_diff = 0 if minelectric==1
					summ m_t_diff if firm_num==`num'
					global firm_t_change = r(mean)
					qui replace t_diff = m_t_diff if firm_num==`num'
					drop m_t_diff

* dummy if technology appears to be cycling
					qui replace T_cycle = abs((T_`scenario' - T_old)/T_old)>$tolerance & (abs((T_`scenario' - T_old3)/T_old3)<0.1*$tolerance | abs((T_`scenario' - T_old4)/T_old4)<0.1*$tolerance) if firm_num==`num'
					qui replace T_`scenario' = T_old if T_cycle==1 & firm_num==`num'
					
* update cycle variables					
					qui replace T_old = T_`scenario' if firm_num==`num'	
					qui replace T_old4 = T_old3 if firm_num==`num'
					qui replace T_old3 = T_old2 if firm_num==`num'
					qui replace T_old2 = T_old if firm_num==`num'					
					
* end fuel consumption rate loop if number of iterations gets large	or cycle variable is on			
					summ num_iter_t if firm_num==`num'
					global num_iter_t = r(mean)
					if $num_iter_t==$max_iter {
						global firm_t_change = 0
					}
					summ T_cycle if firm_num==`num'
					global cycle = r(mean)
					if $cycle>0 {
						global firm_t_change = 0
					}
					
* close fuel consumption rate loop	
					drop firm_wt firm_fcr firm_fcr_old
				}
					
* reset fuel consumption rate tolerance, number of iterations, and cycle variables
				global firm_t_change = 1
				qui replace num_iter_t = 0
				
************
* Loop 4b: iterate until firm's prices converge		
************					
					
* set tolerance
				while $firm_price_change>$tolerance {						
					
* predict new market share and sales			
					drop pr_mkt_share_`scenario' vpr_sales_`scenario'
					sort demo vehicle
					by demo: egen pr_mkt_share_`scenario' = sum(exp(util_new_`scenario'))
					qui replace util_used_`scenario' = ln(used_sh*pr_mkt_share_`scenario'/(1 - used_sh))
					qui replace pr_mkt_share_`scenario' = exp(util_new_`scenario')/(exp(util_used_`scenario') + pr_mkt_share_`scenario')
					qui replace pr_mkt_share_`scenario' = 0 if pr_mkt_share_`scenario'==.
					qui replace pr_sales_`scenario' = pr_mkt_share_`scenario'*mkt_size
					sort vehicle demo
					by vehicle: egen vpr_sales_`scenario' = sum(pr_sales_`scenario')
					
* compute own derivative term
					gen own_deriv = alpha*pr_mkt_share_`scenario'*(1 - pr_mkt_share_`scenario')

* compute cross derivative term
					gen cross_deriv = -alpha*pr_mkt_share_`scenario'
	
* create vectors of own and cross price derivatives
					sort firm_num vehicle demo
					mkmat own_deriv if firm_num==`num'
					mkmat cross_deriv if firm_num==`num'	

* create vector for marginal costs and feebate	
					mkmat mc_`scenario' if firm_num==`num' & demo==1
					mkmat feebate_`scenario' if firm_num==`num' & demo==1		
					mkmat zev_subsidy if firm_num==`num' & demo==1		
					
* create matrix with market shares	
					sort firm_num vehnum demo
					mkmat pr_mkt_share_`scenario' if firm_num==`num'
					matrix define S = J($numveh, $maxdemo, 0)
					forvalues j = 1/$numveh {
						forvalues k = 1/$maxdemo {
							matrix define S[`j',`k'] = pr_mkt_share_`scenario'[`k' + (`j' - 1)*$maxdemo,1]
						}
					}			

* create vector for market size
					mkmat mkt_size if firm_num==`num' & vehnum==1
	
* create matrix of derivatives
					matrix define deriv = J(1,$maxdemo,0)
					matrix define D = J($numveh, $numveh,0)
					forvalues j = 1/$numveh {
						forvalues k = 1/$numveh {
							if `j'==`k' {
								forvalues l = 1/$maxdemo {
									matrix define deriv[1,`l'] = own_deriv[(`j' - 1)*$maxdemo + `l',1]
								}
							}
							if `j'!=`k' {
								forvalues l = 1/$maxdemo {
									matrix define deriv[1,`l'] = cross_deriv[(`k' - 1)*$maxdemo + `l',1]*pr_mkt_share_`scenario'[(`j' - 1)*$maxdemo + `l',1]
								}			
							}
							matrix define D[`j',`k'] = deriv*mkt_size
						}
					}
					
* invert matrix and compute price
					matrix define inv_D = inv(D)
					matrix define p = mc_`scenario' + feebate_`scenario' - zev_subsidy - inv_D*S*mkt_size			

* fill in new prices
					sort firm_num vehnum demo
					forvalues num_m = 1/$numveh {
						qui replace price_new_firm = p[`num_m',1] if firm_num==`num' & vehnum==`num_m'
						qui replace price_new_firm = $price_tol*trans_price  if price_new_firm>$price_tol*trans_price & firm_num==`num' & vehnum==`num_m'
						qui replace price_new_firm = trans_price/$price_tol if price_new_firm<trans_price/$price_tol & firm_num==`num' & vehnum==`num_m'
					}

* update utility and market shares						
					drop util_new_`scenario' pr_mkt_share_`scenario' pr_sales_`scenario' vpr_sales_`scenario'
					gen util_new_`scenario' = util_new + alpha_`scenario'*(price_new_firm - subsidy - trans_price) + beta_cpm*(cpm_`scenario' - cpm) + beta_hw_`scenario'*(ln_hw_`scenario' - ln_hw)
					sort demo vehicle
					by demo: egen pr_mkt_share_`scenario' = sum(exp(util_new_`scenario'))
					qui replace util_used_`scenario' = ln(used_sh*pr_mkt_share_`scenario'/(1 - used_sh))						
					qui replace pr_mkt_share_`scenario' = exp(util_new_`scenario')/(exp(util_used_`scenario') + pr_mkt_share_`scenario')
					qui replace pr_mkt_share_`scenario' = 0 if pr_mkt_share_`scenario'==.
					gen pr_sales_`scenario' = pr_mkt_share_`scenario'*mkt_size
					sort vehicle demo
					by vehicle: egen vpr_sales_`scenario' = sum(pr_sales_`scenario')						
					
* compute price change
					egen m_p_diff = max(abs((price_new_firm - price_old_firm)/price_old_firm)) if firm_num==`num'
					summ m_p_diff if firm_num==`num'
					global firm_price_change = r(mean)
					replace price_old_firm = price_new_firm if firm_num==`num'
					drop own_deriv cross_deriv m_p_diff 					
				}
				
* reset price change
				global firm_price_change = 1				
	
* close loop through firms
			}

* compute difference between fuel consumption rate and standard
			egen wt = sum(pr_sales_`scenario'*multiplier_`scenario') 
			qui replace wt = pr_sales_`scenario'*multiplier_`scenario'/wt
			egen mkt_fcr = sum(wt*fcr_adj_`scenario')
			egen mkt_standard = sum(wt/standard_`scenario')
			gen diff_f = abs((mkt_fcr - mkt_standard)/mkt_standard)
			summ mkt_fcr mkt_standard
			summ diff_f
			global cafe = r(mean)	
	
* adjust shadow price to make up the distance to the standard
			qui replace cafe_price_`scenario' = 1.5*(2 - num_iter_f/$max_iter)*cafe_price_`scenario'*((mkt_fcr - mkt_standard)/mkt_fcr) + cafe_price_`scenario'
			drop wt mkt_fcr mkt_standard diff_f

* end CAFE loop if shadow price gets too high or too low
			gen end = (cafe_price_`scenario'>$cafe_price_change*cafe_price | cafe_price_`scenario'<cafe_price/$cafe_price_change) & num_iter_f>3
			summ end
			global end = r(mean)
			if $end==1 {
			    global cafe = 0				
			}
			drop end
			
* set shadow price to boundary if exceeds boundaries			
			qui replace cafe_price_`scenario' = $cafe_price_change*cafe_price if cafe_price_`scenario'>$cafe_price_change*cafe_price
			qui replace cafe_price_`scenario' = cafe_price/$cafe_price_change if cafe_price_`scenario'<cafe_price/$cafe_price_change

* end CAFE loop at maximum number of iterations
			summ num_iter_f
			global num_iter_f = r(mean)
			if $num_iter_f==0.5*$max_iter {
				global cafe = 0
			}
			qui replace num_iter_f = num_iter_f + 1	
			
* update feebate and markup
			qui replace feebate_`scenario' = cafe_price_`scenario'*(fcr_adj_`scenario' - 1/standard_`scenario')
			qui replace markup_`scenario' = price_new_firm - mc_`scenario' - feebate_`scenario' + zev_subsidy

* close loop on CAFE standard			
		}
		
* reset distance to standard and number of iterations
		global cafe = 1
		qui replace num_iter_f = 0
		
* compute 95th percentile of price changes for entire market, excluding vehicles that were cycling
		qui replace price_new = price_new_firm
		egen m_mkt = max(abs((price_new - price_old)/price_old)) if T_cycle==0 & t_diff>=0 & t_diff<=$tolerance
		egen p95_mkt = pctile(abs((price_new - price_old)/price_old)) if T_cycle==0 & t_diff>=0 & t_diff<=$tolerance, p(95)
		summ m_mkt p95_mkt
		summ p95_mkt
		global price_change = r(mean)
		qui replace price_old = price_new
		drop m_mkt p95_mkt
		
* close loop on market-level price changes		
	}

* reset market-level changes
	global price_change = 1
	
* save simulation outcomes as new variables
	gen sim_price_`scenario' = price_new
	gen sim_used_price_`scenario' = util_used_`scenario'/alpha
	gen sim_T_`scenario' = T_`scenario'
	gen sim_ln_hw_`scenario' = ln_hw_`scenario'
	gen sim_mpg_`scenario' = mpg_`scenario'
	gen sim_fcr_`scenario' = fcr_`scenario'
	gen sim_fcr_adj_`scenario' = fcr_adj_`scenario'
	gen sim_mc_`scenario' = mc_`scenario'
	gen sim_cpm_`scenario' = cpm_`scenario'
	gen sim_cafe_price_`scenario' = cafe_price_`scenario'
	gen sim_feebate_`scenario' = feebate_`scenario'
	drop T_`scenario' ln_hw_`scenario' mpg_`scenario' mc_`scenario' cpm_`scenario' fcr_`scenario' fcr_adj_`scenario' cafe_price_`scenario' feebate_`scenario'
	save "L:\Project-MaritzCX\Workspace1\Public\CAFE Ex Post Replication\Replication Results\table7_simulations", replace
}

*****


log close
